## R code to run models, get tables, and plot figures

## 1. More data preparation ----------------------------------------------------
dat <- dat %>%
  mutate(
    aid_crsc_lag1=ifelse(year>1962&is.na(aid_crsc_lag1), 0, aid_crsc_lag1),
    aid_crsio_lag1=ifelse(is.na(aid_crsio_lag1), 0, aid_crsio_lag1),
    totalaid_lag1=aid_crsc_lag1+aid_crsio_lag1,
    ivvar1 = log((ross_oil_value_2000_lag1 + 1) * (neighbormean_punish_percent_cumsum + 1)+1),
    ivvar2 = log((totalaid_lag1) * (neighbormean_punish_percent_cumsum +1)+1)
  )
  
## 2. Main models -------------------------------------------------------------
## Fixed effect models
fe1 <- lm(LJI ~ threatS + country_name, dat); summary(fe1)
fe2 <- lm(LJI ~ threatS  + gdpS + demdurS + country_name,dat); summary(fe2)
## IV models
iv_total1 <- ivreg(LJI ~ threatS + gdpS + demdurS | gdpS +  demdurS + ivvar1,
                     data=dat); summary(iv_total1)  
iv_total2 <- ivreg(LJI ~ threatS + gdpS + demdurS | gdpS +  demdurS  +  ivvar1 + ivvar2,
                     data=dat); summary(iv_total2)    

## Robust standard errors for fixed effects models ---------------------------
G <- length(unique(dat$country_name))
N <- nobs(fe1)
dfa <- (G/(G-1))* (N-1)/fe1$df.residual
clustered_vcov <- dfa * vcovHC(fe1, type="HC0", cluster="group", adjust=TRUE)
G <- length(unique(dat$country_name))
N <- nobs(fe2)
dfa <- (G/(G-1))* (N-1)/fe2$df.residual
clustered_vcov2 <- dfa * vcovHC(fe2, type="HC0", cluster="group", adjust=TRUE)
fe1.robust.se <- sqrt(diag(clustered_vcov))
fe2.robust.se <- sqrt(diag(clustered_vcov2))
fe1.coeftest <- coeftest(fe1, vcov=clustered_vcov)
fe1.pvals <- fe1.coeftest[,4]
fe2.coeftest <- coeftest(fe2, vcov=clustered_vcov2)
fe2.pvals <- fe2.coeftest[,4]

## Robust se for IV models -----------------------------------------------------
# texreg does not by default include IV gof statistics
G <- length(unique(dat$country_name))
N <- nobs(iv_total1)
dfa <- (G/(G-1))* (N-1)/iv_total1$df.residual
clustered_vcov_ivtotal1 <- dfa * vcovHC(iv_total1, type="HC0", cluster="group", 
                                        adjust=TRUE)
V1.robust.se <- sqrt(diag(clustered_vcov_ivtotal1))
V1.coeftest <- coeftest(iv_total1, vcov= clustered_vcov_ivtotal1)
V1.pvals <- V1.coeftest[,4]
tr <- texreg::extract(iv_total1,include.rsquared=FALSE,include.adjrs=FALSE)
tr@gof <- c(tr@gof, summary(iv_total1,diagnostics=TRUE)$diagnostics[1,4])
tr@gof.names <- c(tr@gof.names, "Weak instruments test")
tr@gof.decimal <- c(tr@gof.decimal, TRUE)
tr@gof <- c(tr@gof, summary(iv_total1,diagnostics=TRUE)$diagnostics[2,4])
tr@gof.names <- c(tr@gof.names, "Wu-Hausman test")
tr@gof.decimal <- c(tr@gof.decimal, TRUE)
tr@gof <- c(tr@gof, summary(iv_total1,diagnostics=TRUE)$diagnostics[3,4])
tr@gof.names <- c(tr@gof.names, "Sargan test")
tr@gof.decimal <- c(tr@gof.decimal, TRUE)
## iv total 2 
N <- nobs(iv_total2)
dfa2 <- (G/(G-1))* (N-1)/iv_total2$df.residual
clustered_vcov_ivtotal2 <- dfa2 * vcovHC(iv_total2, type="HC0", cluster="group", 
                                         adjust=TRUE)
V2.robust.se <- sqrt(diag(clustered_vcov_ivtotal2))
V2.coeftest <- coeftest(iv_total2, vcov= clustered_vcov_ivtotal2)
V2.pvals <- V2.coeftest[,4]
tr2 <- texreg::extract(iv_total2,include.rsquared=FALSE,include.adjrs=FALSE)
tr2@gof <- c(tr2@gof, summary(iv_total2,diagnostics=TRUE)$diagnostics[1,4])
tr2@gof.names <- c(tr2@gof.names, "Weak instruments test")
tr2@gof.decimal <- c(tr2@gof.decimal, TRUE)
tr2@gof <- c(tr2@gof, summary(iv_total2,diagnostics=TRUE)$diagnostics[2,4])
tr2@gof.names <- c(tr2@gof.names, "Wu-Hausman test")
tr2@gof.decimal <- c(tr2@gof.decimal, TRUE)
tr2@gof <- c(tr2@gof, summary(iv_total2,diagnostics=TRUE)$diagnostics[3,4])
tr2@gof.names <- c(tr2@gof.names, "Sargan test")
tr2@gof.decimal <- c(tr2@gof.decimal, TRUE)
## save the table  ------------------------------------------------------------
texreg(list(fe1, fe2, tr, tr2),
       file = dir_table,
       override.se = list(fe1.robust.se, fe2.robust.se, V1.robust.se, V2.robust.se),
       override.pvalues= list(fe1.pvals, fe2.pvals, V1.pvals, V2.pvals),
       caption="\\textbf{De facto independence and the demand for insurance.} 
       Four models of de facto independence and the demand for insurance. 
       Models 1 and 2 are linear models. Model 3 and 4 are 2SLS models instrumenting 
       for the demand for insurance (Model 4 with multiple instruments). All models
       include robust standard errors clustered by country and country fixed effects 
       (not shown).",
       omit.coef="name",caption.above=TRUE,
       booktabs=TRUE,dcolumn=TRUE,use.packages=FALSE,
       label="table:iv1",
       custom.coef.names=c("Intercept","Demand for insurance","log(GDP/capita)",
                           "Years democratic (logged)"),
       include.aic=FALSE,include.deviance=FALSE,
       include.rmse=FALSE,style="apsr",digits=2,scalebox=0.75, float.pos="t!") 

## Plot the marginal effects  -------------------------------------------------
iv_plot <- ivreg(LJI ~ threat +  gdpS + demdurS  | gdpS +  demdurS  + ivvar1 + ivvar2,
                 data=dat)    
dfaP <- (G/(G-1))* (N-1)/iv_plot$df.residual
clustered_vcov_ivplot <- dfaP * vcovHC(iv_plot, type="HC0", cluster="group", 
                                       adjust=TRUE)
beta <- c(coef(iv_plot)[1:4])
covvar <- clustered_vcov_ivplot[1:4,1:4]
x <- cbind(1,
           seq(min(dat$threat, na.rm = TRUE),max(dat$threat, na.rm = TRUE),length.out=101), 
           mean(dat$gdpS, na.rm=TRUE),
           mean(dat$demdurS, na.rm=TRUE))
beta.sim <- mvrnorm(10000,beta,covvar)
xb <- x %*% t(beta.sim) %>% pmax(., 0) %>% pmin(., 1) ## Note here that LJI is bounded between 0 and 1
confinf <- apply(xb,1,quantile,probs = c(0.025, 0.975)) 
plot_dat <- data.frame(threat = x[,2], 
                       lower = confinf[1,], 
                       middle = apply(xb,1,mean,na.rm=TRUE), 
                       upper = confinf[2,])
plot_me <- ggplot(plot_dat, aes(x = threat)) +
  geom_ribbon(aes(ymin = lower, ymax= upper),
              fill="gray70") + 
  scale_x_continuous('Demand for insurance',
                     breaks=c(0,0.2,0.4,0.6,0.8,1)) + 
  scale_y_continuous('E(Judicial independence)' )	+ 
  geom_line(aes(y = middle), linetype = 1,size=1.3) + 
  theme_bw() + 
  theme(axis.title.y = element_text(size=14,angle=90),
        axis.text.y=element_text(size=14)) + 
  theme(axis.title.x = element_text(size=14),
        axis.text.x=element_text(size=14))